ELSEVIER 


Available  online  at  www.sciencedirect.com 

ScienceDirect 

Journal  of  Power  Sources  178  (2008)  269-281 

www.elsevier.com /locate /jpowsour 


JOURNAL  OF 


Three-dimensional,  single-phase,  non-isothermal 
CFD  model  of  a  PEM  fuel  celL 

Carlos  Martinez  Bacaa *  *,  Rowland  Travis a,  Mads  Bangb 

a  Department  of  Mechanical  Engineering,  Imperial  College,  London  Exhibition  Road,  South  Kensington,  London  SW7  2AZ,  UK 
b  Institute  of  Energy  and  Technology,  Aalborg  University,  Aalborg,  Denmark 

Received  29  August  2007;  received  in  revised  form  14  November  2007;  accepted  1  December  2007 


Abstract 

A  comprehensive,  three-dimensional  analysis  of  a  polymer  electrolyte  membrane  (PEM)  fuel  cell  has  been  developed  to  study  the  performance  of 
this  device  under  different  operational  conditions.  This  steady-state  analysis  is  single-phase  and  non-isothermal.  A  commercial  computational  fluid 
dynamics  (CFD)  program  provided  a  numerical  platform  for  solving  the  conservation  equations  for  species,  energy,  charge,  mass  and  momentum. 
Different  boundary  conditions  were  added  to  a  computational  domain  to  simulate  single  channel  PEM  fuel  cell.  The  electrochemistry  involved  in 
this  model  was  added  by  a  set  of  user-defined  subroutines  that  feature:  electrochemical  reactions,  electric  and  ionic  charge  and  heat  generation. 
The  calculations  were  then  solved  by  an  iterative  method  following  an  adapted  computational  procedure.  The  results  were  validated  with  other 
computational  models  and  experimental  data.  These  show  a  noticeable  non-uniform  distribution  of  the  current  density  across  the  catalyst  layer 
(CL)  at  different  operational  conditions.  The  results  emphasize  on  the  differences  of  anodic  and  cathodic  activation  overpotentials,  the  oxygen 
transport  limitations  and  the  ohmic  losses  distributions  of  both  proton  and  electric  overpotentials. 
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1.  Introduction 

A  fuel  cell  is  as  an  electrochemical  converter.  It  can  be  said 
that  a  fuel  cell  is  a  device  that  converts  chemical  energy  into 
electrical  power  and  ideally  can  continue  to  operate  as  long  as  it 
is  fed  with  suitable  fuels  and  oxidants  and  the  reaction  products 
are  being  removed. 

Electrochemical  power  sources  differ  from  others,  such  as 
thermal  power  plants  because  the  energy  conversion  occurs  with¬ 
out  any  intermediate  steps.  For  example,  in  the  case  of  thermal 
power  plants,  fuel  is  first  converted  into  thermal  energy  and 
then  into  electric  power  via  generators.  In  fuel  cells  systems  this 
multi-step  process  is  achieved  by  electrochemical  reactions.  As 
a  consequence,  electrochemical  systems  show  some  advantages, 
such  as  energy  efficiency. 
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There  is  extensive  literature  available  on  numerical  ways  of 
modelling  fuel  cells.  Polymer  electrolyte  membrane  (PEM)  fuel 
cells  have  been  modelled  for  more  than  15  years.  Although  one¬ 
dimensional  in  nature  some  of  the  early  analyses  remain  the 
basis  for  most  of  the  more  elaborated  models  to  date.  Bernardi 
and  Verbrugge  [5]  made  a  big  emphasis  on  the  electrochemistry 
and  diffusive  processes  inside  the  fuel  cell’s  membrane  electrode 
assembly  (ME A).  Springer  et  al.  [22]  provided  a  great  insight 
into  the  water  transport  process  inside  the  MEA.  Ticianelli  et 
al.  [25]  made  a  numerical  model  which  introduced  a  way  of 
calculating  locally  the  cell  potential.  Fuller  and  Newman  [11] 
produced  a  model  that  could  predict  the  performance  of  the  cell 
along  a  longitudinal  channel  by  integrating  the  solution  at  var¬ 
ious  points  down  the  channel.  They  pinpointed  the  effect  of 
dehydration  in  the  fuel  cell’s  membrane  due  to  the  electroos¬ 
motic  drag.  Lampimen  and  Fomino  [15]  introduced  entropy 
changes  in  to  their  model.  Gurau  et  al.  [12]  used  a  special  han¬ 
dling  of  the  transport  equations  which  enabled  them  to  use  the 
same  numerical  method  for  a  unified  computational  domain. 
This  domain  constituted  a  polymer  membrane,  two  catalyst  lay¬ 
ers  (CLs),  two  gas  diffusers  (GDs)  and  two  feeding  channels. 
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Fig.  1.  Diagram  of  the  solution  procedure. 


This  in  turn  eliminated  the  need  for  prescribing  arbitrary  bound¬ 
ary  conditions  for  the  fuel  cell  interfaces.  The  papers  mentioned 
above  are  the  foundation  of  the  model  presented  in  this  paper. 

As  PEM  fuel  cell  models  were  becoming  more  specific,  peo¬ 
ple  began  to  realize  that  one  of  the  most  complicated  issues  in 
their  analyses  was  the  transport  of  water  in  the  PEM.  Zawodzin- 
ski  et  al.  [28]  experimented  on  several  types  of  proton  exchange 
membranes.  They  reported  water  sorption  characteristics,  dif¬ 
fusion  coefficient,  electroosmotic  drag  and  conductivity.  Van 
Zee  et  al.  [16]  and  Shimpalee  et  al.  [17]  conducted  a  thorough 
experiment  concerning  water  transport  in  the  membrane  that  was 
validated  with  a  numerical  model.  Most  of  their  papers  conclude 
that  the  velocity  of  the  water  inside  the  membrane  is  too  small 
and  has  little  effect  on  the  overall  water  mass  flow  rate.  Later 
on  Zawodzinski  et  al.  [27]  and  more  recently  Kulikovsky  [14] 
and  Sui  and  Djilali  [23]  confirmed  that  the  electroosmotic  drag 
coefficient  can  be  considered  as  constant  for  a  wide  range  of 
water  content  in  the  membrane. 

The  three-dimensional  PEM  fuel  cell  model,  presented  by 
Dutta  et  al.  [9]  showed  how  to  modify  a  commercial  compu¬ 
tational  fluid  dynamics  (CFD)  code  to  include  the  necessary 
electrochemical  processes.  This  and  subsequent  publications 
[10,18]  stressed  the  importance  of  the  diffusive  and  convective 
transport  of  reactants  and  products  in  the  MEA  and  feeding  chan¬ 
nels.  These  publications  also  introduced  geometric  changes  in 
their  analyses.  This  approach  derived  in  a  better  understanding  of 
the  overall  geometry  of  the  cell  and  gave  way  to  numerous  design 
analyses.  Later,  some  papers  focused  exclusively  on  how  to  opti¬ 
mise  the  dimensions  of  specific  PEM  fuel  cell  components  such 
as  gas  channel  aspect  ratios,  gas  diffusers  and  electrolyte  thick¬ 
nesses,  nation  content  on  CLs  and  catalyst  porosities  (Fig.  1). 

Most  PEM  fuel  cell  models  available  in  the  literature  assume 
that  the  electrochemical  reactions  take  place  at  an  interface 
between  the  diffusers  and  the  electrolyte.  It  is  convenient  to 
model  CLs  as  interfaces  because  it  is  possible  to  define  arbitrary 
boundary  conditions  at  both  surfaces  of  this  “wall”  such  as  heat 
transfer,  mass  sources  or  fixed  temperature  and  concentrations. 
In  reality  CLs  are  normally  a  few  microns  thick  and  in  this  paper 
we  address  the  distributions  of  the  activation  overpotentials  and 
current  densities  across  the  CLs.  The  computational  procedure 
presented  in  this  model  is  able  to  measure  spatial  distribution 
of  these  losses  inside  these  thin  layers.  In  general  we  try  to  pre¬ 
dict  the  distributions  of  the  overpotentials  across  the  model’s 
MEA. 

2.  Model  description 

2.7.  Fuel  cell  principles 

The  computational  domain  of  our  model  is  shown  in  Fig.  2. 
The  model  constitutes  the  anode,  cathode  and  electrolyte  regions 
of  a  PEM  fuel  cell.  Hydrogen-rich  fuel  is  fed  in  to  the  anode 
channel,  where  part  of  this  travels  across  the  GD.  At  the  anode 
CL  an  oxidation  reaction  splits  the  hydrogen  into  protons  and 
electrons  according  to: 

H2  -*  2H+  +  2e~ 


(1) 
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Fig.  2.  Computational  domain  of  the  model. 


Driven  by  an  electric  field,  the  protons  migrate  across  the 
PEM  from  the  anode  towards  the  cathode.  The  free  electrons 
are  conducted  across  the  cell  normally  through  the  platinum 
supported  carbon  black  particles  embedded  in  the  CLs,  then 
through  the  electric  conductive  fibers  of  the  GDs,  the  current 
collectors  (CCs)  and  through  an  external  circuit  which  may 
include  a  load  and  back  again  finishing  at  the  cathode  CL  (see 
Scheme  1). 

The  cathode  channels  are  fed  with  air.  The  oxygen  in  the  air 
travels  through  the  GD  towards  the  CL  where  it  combines  with 
the  migrating  protons  and  the  electrons  to  form  water  according 
to  the  following  exothermic  reduction  reaction: 

^02  +  2H+  +  2e~  ->  H20  (2) 

The  overall  process,  consumes  reactants  generating  electric¬ 
ity,  water  and  heat. 

The  steady- state  PEM  fuel  cell  model  presented  in  this  paper 
is  three-dimensional,  single-phase  and  non-isothermal  and  it 
comprises  the  following  features: 


“electroneutral”  and  behaves  as  a  gas  barrier.  It  is  also  assumed 
that  there  is  a  constant  potential  at  the  interface  between  the 
anode  GD  and  the  CCs.  Lurther,  these  interfaces  and  the  chan- 
nels/CCs  interfaces  are  assumed  to  be  adiabatic.  We  think  that 
the  biggest  limitation  of  this  model  is  the  assumption  of  single¬ 
phase  modelling. 

2.2.  Model's  domain  and  geometry 

Lig.  2  shows  the  model’s  computational  domain.  It  contains 
around  30,000  nodes  divided  into  four  separate  domains.  These 
are:  channels,  GDs,  CLs  and  a  PEM. 

The  channels  are  12  mm  long  and  have  a  cross  sectional 
area  of  0.25  mm2.  The  MEA  represents  the  GDs,  CLs  and  elec¬ 
trolyte.  It  is  10  mm  long,  1  mm  wide  and  0.69  mm  thick.  The 
GD,  CL  and  electrolyte  are  300,  20  and  50  pirn  thick,  respec¬ 
tively. 


2.3.  Electrochemical  reaction  rate 

The  local  current  densities,  /a  and  zc,  can  be  accurately 
predicted  at  each  of  the  catalyst  layers  elements  according 
to  Bard  and  Laulkner  [3].  The  correct  implementation  of  the 
Butler- Volmer  equation  in  the  model  gives  the  local  current 
density  at  the  cathode  (Eq.  (3))  and  anode  (Eq.  (4))  catalyst 
layers. 

4(l-ac)F^ 

rj,  ^act,c 

(3) 


lc  —  io,cXo2 


exp 


4  acF 


RT 


'  *7act,c 


•  multi-component  flow, 

•  convective  and  conductive  heat  transfer, 

•  transport  of  species  across  a  porous  media, 

•  electrochemical  reactions, 

•  electric  and  protonic  charge  conduction. 

The  equations  governing  these  processes  include  the  full 
mass,  momentum,  species,  charge  and  energy  conservation 
equations.  Given  some  appropriate  boundary  conditions,  an 
extensive  suite  of  user  subroutines  and  a  costumed  itera¬ 
tive  procedure,  the  commercial  CFD  solver  Star-CD  v3.24 
provided  the  numerical  platform  to  process  the  calcula¬ 
tions. 

This  model  assumes  that  the  gases  in  the  channels,  GDs  and 
CLs  behave  as  compressible  ideal  gases  and  that  the  membrane  is 


G  —  ^0,a^H2 


exp 


exp 


(  2(1  —aa)F  V 
V  RT  llacUa)_ 
(4) 


where  /o,c/a>  ^act,c/a  and  aa/c  are  the  exchange  current  densi¬ 
ties,  the  activation  overpotentials  and  the  transfer  coefficients, 
respectively.  Subscripts  a  and  c  refer  to  anode  and  cathode, 
respectively.  The  two  unknown  variables  in  these  equations 
are  the  concentration  of  the  reactants  and  the  activation  over¬ 
potentials.  The  concentration  of  the  reactants  is  solved  by 
the  conservation  of  species  equation  in  the  catalyst  lay¬ 
ers  (see  Section  2.4.3).  The  activation  overpotentials  have  a 
special  solving  procedure  which  is  given  below  in  Section 
2.5. 


2.4.  Model  equations 


H20 1  Ib2  k*e- 

t  2T 

I  H2<H 

_ _ I 


cathode  channel 
^athode  GDL 
ithode  CL 
electrolyte 
anode  CL 
^•anode  GDL 
anode  channel 


This  model  has  been  programmed  to  work  with  two  gas  mix¬ 
tures,  hydrated  hydrogen  and  hydrated  air.  The  transport  of  gas 
mixtures  is  solved  by  the  mass  (Eq.  (5)),  momentum  (Eq.  (6)) 
and  energy  (Eq.  (11))  conservation  equations: 

2-^)  =  *, 


Scheme  1 .  Schematic  PEM  fuel  cell. 


(5) 
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d 

—  {pUjUi 


(6) 


where  sp  is  the  mass  source  term,  x/  are  the  cartesian  coordinates 
(i  =  x,  y,  z),Ui  the  absolute  fluid  velocity  component  in  direc¬ 
tion  Xi,  x ij  the  viscous  stress  tensor  components,  P  the  pressure 
and  p  the  density. 

The  viscous  stress  tensor  component  r/y  in  Eq.  (6)  is  calcu¬ 
lated  as  follows: 


Tij  =  [l 


f  dui  duj\ 
Vdxy  +  dxi  ) 


duk 

- 6ii 

dxk 


(7) 


where  /x  is  the  molecular  viscosity,  Sq  the  Kronecker  delta. 

Each  species  m  of  the  gas  mixtures,  whose  local  concentration 
is  expressed  as  a  mass  fraction  Ym,  is  assumed  to  be  governed 
by  a  species  conservation  equation  of  the  form: 


d 

dXj 


OmP 


dTm 

dXj 


Sm 


(8) 


where  sm  the  rate  of  mass  production  of  species  m  and  Dm  is 
the  molecular  diffusivity  of  species  m. 

The  molecular  diffusivity  of  the  gas  mixtures  changes  due  to 
temperature,  pressure  and  species  concentrations.  During  oper¬ 
ation,  fuel  cells  vary  their  gas  compositions  due  to  chemical 
reactions.  Therefore  a  variable  diffusion  coefficient  for  the  gases 
is  needed. 

In  a  binary  gas  mixture  with  components  m  and  n ,  such  as 
hydrogen  and  water  vapour,  the  binary  diffusion  coefficients, 
Dmn ,  can  be  approximated  by  Eq.  (9)  according  to  Cussler  [8]: 


D 


mn 


(9) 


where  Dmn^ f  is  the  reference  diffusion  coefficient  of  the  mix¬ 
ture  at  pressure  Pre f  and  temperature  Tref .  In  ternary  systems,  the 
calculation  of  effective  diffusion  coefficients  can  become  very 
complex.  Nevertheless  PEM  fuel  cells  operate  at  a  range  of  tem¬ 
peratures  and  pressures  which  allow  alternative  approximation 
methods  to  these  coefficients  to  be  as  accurate  as  the  traditional 
and  more  elaborate  models.  Mills  [20]  suggested  a  relation  based 
on  a  “mixture  rule”  that  proves  applicable  to  PEM  fuel  cell  gas 
mixtures  [2],  such  as  oxygen,  nitrogen  and  water  vapour,  and  is 
given  as 

1  -  Xm 

An, mix  =  - - -  (10) 

J2  ( Xn/Dmn ) 

n=l  \nfm 


equation. 


9  /  ,  ,8r 

7j—  |  ph\U j  k  —  ^  /  jht,mPVm,j 


dXj 


d  P  dui  v-^ 

=  Ujdx-+  T,; aTT  +  Sh  -  2^Hn 


dXi 


(11) 


Here  sm  is  the  rate  of  production  or  consumption  of  species 
m  due  to  a  chemical  reaction,  k  is  the  thermal  conductivity, 
Vm,j  the  diffusion  velocity  vector  in  the  j  direction  and 
the  thermal  enthalpy  of  species  m  which  is  defined  by  Eq. 
(11) 


—  CpT  cpjQfTref  (12) 

where  cp  is  the  mean  constant-pressure  specific  heat  at  tem¬ 
perature  T  and  cAref  a  reference  specific  heat  at  temperature 
Tref- 

The  reactants  and  products  are  assumed  to  behave  as  ideal 
gases,  thereby  the  density  p  of  the  gas  mixture  is  calculated 
with 

P=  _  -  (13) 

m 

where  Mm  is  the  molecular  weight  of  species  m. 


2.4.1.  Conservation  equations  in  the  channels 

The  transport  of  gaseous  species  in  the  fuel  cell  channels,  is 
modelled  with  the  following  conservation  equations  for  mass 
(Eq.  (14)),  momentum  (Eq.  (15)),  species  (Eq.  (16))  and  energy 
(Eq.  (17)) 


s; {puj) = 0 

3  ,  ,  dp 

—  (PUjUi  -  X ij)  =-  — 


dx 


3 Xi 


(14) 


(15) 


9  /  dYm 

|  pu j Ym  —  Dm  m\xp^—  J  —  0 


dx 


(16) 


3  /  ,  ,  dT  ' 

dx  '  P^tUj  Xmix  ^  +  y  hj,m pVm,j 


‘dX; 


dp 
dx  ■ 


dUr 

+  Tijd7; 


(17) 


where  Xm  refers  to  the  molar  fraction  of  species  m.  Eq.  (10) 
calculates  the  effective  diffusion  coefficient,  Dm?m ix  of  species 
m  in  a  mixture  composed  of  several  different  components.  For 
further  insight  into  the  validation  and  background  of  Mills’s 
model  please  refer  to  [1]. 

Heat  transfer  is  implemented  with  Eq.  (11).  This  is  obtained 
by  multiplying  the  species  conservation  equation  (Eq.  (8)) 
of  each  species  m  of  the  mixture  by  its  heat  of  forma¬ 
tion,  Hm,  and  subtracting  it  from  the  energy  conservation 


where  km\x  is  the  mixture  heat  conductivity  of  species  m  given 
by  Eq.  (18)[19]  formula 


.eff 

mix 


1 

2 


(18) 


where  kn  is  the  thermal  conductivity  of  component  n. 

The  velocities  at  the  inlet  boundary  of  the  channels  were 
set  as  a  function  of  the  total  current  produced,  /,  the  channel 
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cross  sectional  area,  Ac h  and  the  stoichiometric  ratio  £.  The 
stoichiometry  ratio,  f ,  is  defined  as 


mass  of  fuel  or  oxidant  input  to  the  cell 

f  = - - -  (19) 

mass  of  fuel  or  oxidant  reacted  in  the  cell 

The  total  current  is  a  variable  dependent  on  the  cell  overpoten¬ 
tials  and  species  concentrations  that  change  during  the  solution 
procedure.  Therefore  velocities  are  updated  on  an  iterative  basis, 
so  that  the  stoichiometric  ratios  at  the  inlet  boundaries  are  always 
kept  constant.  These  are  given  as 


uc  —  £c 


Mq2I 


Yn ,  in  Or  in  Arh 


Ma-(a, 


Mu2I 


(20) 


where  uc  and  wa  are  the  inlet  velocities  and  £c  and  £a  the  stoi¬ 
chiometry  ratios  for  cathode  and  anode,  respectively.  Mo2  and 
Mh2  are  the  molar  masses  for  oxygen  and  hydrogen.  The  mass 
fractions  of  oxygen  and  hydrogen  at  the  inlets  are  represented 
hy  To2,in  and  Tj^in  and  the  densities  for  both  gas  mixtures,  at 
the  cathode  and  anode  inlet  boundaries  are  given  by  pc?in  and 
pa?in,  respectively.  At  the  channel  outlet  boundaries,  Neumann 
boundary  conditions  (zero  concentration  gradient  normal  to  the 
outlet  face)  are  prescribed  for  the  normal  velocities,  tempera¬ 
ture  and  species  equations.  Pressure  is  prescribed  at  the  outlet 
boundaries  to  match  the  desired  operational  pressure.  The  tem¬ 
perature  at  the  surrounding  walls  of  the  channels  is  assumed  to 
be  constant. 


2.4.2.  Conservation  equations  in  the  gas  diffusers 

The  gas  transport  in  the  GDs  is  restricted  by  the  material’s 
porosity,  sq d,  which  is  the  fractional  void  volume.  This  changes 
the  conservation  equations  in  the  following  way: 


—  (pSGDUj)  =  0 


dx 


u  j  =  —- 


kp,GD  dp 


dx ; 


fl 


dx  i 


—  I  pSGDUjYm  -  DefGDpsGD—^ 


dY„ 


dXi 


(21) 

(22) 

(23) 


9  (  dT  ^ 

—  [phteGDU  j  -  &gd(1  -  £gd)—  +  'ffKmpVmj 


dP  dUi 

=  M  7 - h  Tij - 

J  dXj  J  dXj 


(24) 


The  mass,  momentum,  species  and  energy  equations  (Eqs. 
(21)-24)  are  reduced  by  £gd-  The  momentum  equation  has 
been  replaced  by  Darcy’s  law,  where  kp$ d  is  the  hydraulic 
permeability  of  the  gas  mixture  in  the  GDs.  The  species  binary 
diffusion  coefficients  in  Eq.  (23)  are  reduced  by  the  Bruggeman 
correction  [4] 

Dfn  =  Dmne®§  (25) 

where  y  is  the  tortuosity  of  the  medium. 

Electric  conduction  through  the  GD  fibers  is  modelled  with 
Eq.  (26).  It  relates  the  flux  of  current  to  the  gradient  of  the 


electric  potential  0e- .  In  this  conduction  process  the  charge  is 
accelerated  by  an  electric  field  and  the  intrinsic  resistance  of  the 
conductive  material 


d 

d.xj 


^“^e-,GD 


dXj 


=  0 


(26) 


where  /cc-?gd  is  the  electric  conductivity  of  the  GD. 

It  is  assumed  that  the  temperature  of  the  CCs  is  uniform. 
Thereby  a  fixed  temperature  was  defined  at  these  boundaries. 
It  is  also  assumed  that  the  electric  potential  at  this  interface  is 
fixed  to  a  reference  value,  (see  Section  2.5). 


2.4.3.  Conservation  equations  in  the  catalyst  layers 

In  this  region  there  are  five  transport  processes:  gas 
species  transport  through  the  pores,  water  transport  through  the 
hydrophilic  electrolyte  regions  and  void  paths,  heat  conductiv¬ 
ity  through  the  solid  matter,  proton  conductivity  through  the 
electrolyte  fraction  and  electric  conductivity  through  the  carbon 
black  particles. 

All  electrochemical  reactions  in  a  PEM  fuel  cell  occur  at 
these  thin  layers.  A  pseudo-homogeneous  model  of  the  CLs  has 
been  adopted  from  [6].  It  has  been  assumed  that  all  materials 
and  properties  of  the  CL  have  been  homogeneously  dispersed 
throughout  its  whole  extent.  Therefore,  general  values  can  be 
prescribed  for  the  entire  catalyst  region. 

The  electric  potential  loss  in  the  PEM  is  related  to  the  fact  that 
an  electric  field  is  necessary  in  order  to  maintain  the  motion  of 
the  protons  through  the  membrane.  This  field  is  provided  by  the 
existence  of  a  potential  gradient  across  the  cell,  which  is  directed 
opposite  from  the  outer  field  that  gives  the  cell  potential,  and  thus 
has  to  be  subtracted.  It  can  be  shown  that  this  loss  obeys  Ohm’s 
law,  [21]:  hence,  it  can  also  be  assumed  that  the  transport  of 
protons  (Eq.  (27))  is  driven  by  a  conduction  process  similar  to 
the  electric  conductivity. 

~~Lj{au+eV&d^~)=su+  (27) 

In  Eq.  (27)aH+  is  the  proton  conductivity,  which  is  not 
constant  and  sH+  the  source  term  for  proton  production. 
The  proton  conductivity  is  affected  by  the  CL  porosity,  £cl- 
A  similar  relation  to  the  Bruggerman  correction  is  used  to 
calculate  the  effective  conductivity  for  protons  in  the  electrolyte 
fraction  of  the  CLs  [2]. 

The  source  terms  of  the  proton  conservation  equation  (Eq. 
(27)),  are  defined  as  the  rate  of  production  of  protonic  charge  per 
unit  volume.  This  scalar  is  equal  to  the  divergence  of  the  current 
density,  given  by  the  current  density  divided  by  the  thickness  of 
the  electrode  on  both,  anode  and  cathode  CLs,  shown  here  in 
Eqs.  (28) 


ic 

%+,c  —  ’  £H+,a 

tCL 


l  Si 

tCL 


(28) 


where  tec  is  the  thickness  of  the  catalyst  layer.  The  electric 
conductivity  in  the  CLs  is  modelled  in  the  same  way  as  in 
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the  GDs 

d 

dxj 


~Ke-,CL 


#e~ 

dxj 


—  Sp. 


(29) 


where  kq-  cl  refer  to  the  CLs  electric  conductivity,  which  is 
assumed  to  be  constant.  Similarly,  the  rate  of  production  of  elec¬ 
tric  charge  given  in  Eq.  (29)  at  the  anode  and  cathode  CLs,  follow 
those  of  the  protons 


l  c 

tcc 


Sp  : 


tc  L 


(30) 


The  conservation  equations  for  mass,  momentum,  species 
and  energy  in  the  CLs  are: 


—  (pecLUj)  =  sp 


Uj  =  — 


fcp.CL  dp 
/I  dxj 


dx; 


—  (  PSCLUjYm  +  ~  OfcLP£CL—^ 


dYlr, 


dx ; 


(31) 


(32) 


(33) 


9  (  dT 

7^7  (  Phl£CLU j  -  ^cl(1  -  eCL)^7  +  y'jh.mPvm,j 
1  \  1  m 


dP  dui  ^ — -v 

=  Uj  —  I”  xij~7  I ~  sh  —  / 

dXj  ^ 


'  dxj 


(34) 


The  binary  diffusion  coefficients  for  gases,  D^CL,  in  the  CLs 
are  calculated  in  a  similar  fashion  as  in  the  GD  (Eq.  (23))  using 
£CL  instead.  The  source  terms  for  the  mass  conservation  equation 
in  the  cathode  and  anode  catalyst  layers  (Eq.  (31))  are  given  by 
Eqs.  (35)  and  (36),  respectively 


p,c  — 


{  4^h2Q 
V  2  F 


+  ^H2Oabs,c 


P,  a  — 


+  ^H2Oabs,a 


(35) 

(36) 


where  Mo2,  Mh2  and  Mh2o  are  the  molar  masses  of  oxygen, 
hydrogen  and  water,  respectively.  The  terms  where  the  current 
density  is  divided  by  the  thickness  of  the  catalyst  layer  in  Eqs. 
(35)  and  (36)  refer  to  the  rate  of  production  of  electric  charge 
per  unit  volume. 

The  source  terms  for  the  oxygen  and  hydrogen  conservation 
equation  (Eq.  (33))  in  the  CLs  are  given  Eqs.  (37)  and  (38), 
respectively 


so2  = 


^h2  = 


(37) 

(38) 


The  main  heat  source  in  PEM  fuel  cells  is  given  by  the  reduc¬ 
tion  reaction  at  the  cathode  CL.  The  heat  sources  of  the  CLs 


are  given  by  Eqs.  (39)  and  (40)  for  the  cathode  and  anode, 
respectively. 


sh,  c  — 


n- as) 

AF 


?7 act, c 


(39) 


sh,a  — 


Cf\ 


eff 


H+ 


(40) 


where  is  the  entropy  change  of  the  oxygen  reduction  reaction. 
The  first  term  on  the  right-hand- side  of  Eq.  (39)  represents  the 
energy  loss  due  to  entropy  changes  as  well  as  irreversibilities 
with  charge  transfer  [15].  The  hydrogen  oxidation  reaction  heat 
source  is  very  small  and  has  little  effect  in  the  cell  performance 
[13]. 


2.4.4.  Conservation  equation  for  water  in  the  electrolyte 

The  electrolyte  in  a  PEM  fuel  cell  is  effectively  a  thin  polymer 
membrane.  It  comprises  the  PEM  and  a  fraction  of  the  CLs 
(nation  content  in  the  CLs).  Its  properties  allow  conduction  of 
heat,  water  and  protonic  charge. 

Water  plays  an  important  role  in  the  overall  cell  performance 
because  PEMs  have  to  be  hydrated  in  order  to  improve  their  pro¬ 
tonic  conductivity.  Water  distribution  may  hinder  or  enhance  the 
fuel  cell  performance.  On  the  one  hand  there  is  a  risk  of  having 
excess  water  which  will  obstruct  the  transport  of  reactants  to  the 
catalyst  sites  due  to  the  agglomeration  of  water  in  the  porous 
medium.  On  the  other  hand  de-humidification  of  the  membrane 
may  occur,  reducing  the  protonic  conductivity.  In  extreme  cases 
of  complete  drying  of  the  membrane,  local  burn-out  of  the  elec¬ 
trolyte  may  take  place  and  hydration  strain  becomes  an  adverse 
effect  threatening  the  integrity  of  the  material. 

The  conservation  equation  for  water  in  the  electrolyte  is  given 
as: 

2-7  (/h2oj)  =  0  (41) 

where  /h2o ,7  is  the  diffusional  flux  of  water  inside  the  polymer. 
There  are  two  driving  forces  that  transport  the  water  absorbed 
in  the  PEM,  diffusion  and  electroosmotic  drag. 

The  PEM  absorbs  water  from  the  surrounding  gases  via  the 
hydrophilic  electrolyte  fraction  in  the  CLs.  This  happens  when 
the  equilibrium  water  content  “(A.)”  of  this  fraction  is  smaller 
than  that  of  the  surrounding  gases.  If  on  the  contrary,  the  X  of  the 
electrolyte  fraction  is  higher  than  that  of  the  surrounding  gases 
then  water  will  be  given  from  the  polymer  to  the  surroundings. 
X  is  defined  as  the  ratio  of  the  number  of  water  molecules  to  the 
number  of  charged  HSO3  sites.  Springer  et  al.  [22]  produced  the 
following  relationship  for  X  based  on  experiments  performed  on 
Nation  117  membranes: 

X  =  0.043  +  17.81a  -  39.85a2  +  a3  a  <  1  (42) 

X  =  14  +  1.4(a  —  1)  a  >  1 


(43) 
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where  the  activity  of  water,  a ,  is  related  to  the  partial  pressure 
of  water  and  the  saturation  pressure  of  water  Psat,  according  to: 

«  =  ,44) 

^sat 

where  Xn2o  is  the  molar  fraction  of  water.  The  saturation  pres¬ 
sure  of  water  Psat  is  given  by  expression  (45)  as  reported  by 
Springer  et  al.  [22] 

logio^sat  =  -2.1794  +  0.02953 T  -  9.1837  x  lO^r2 

+  1.4454  x  10“7r3  (45) 


where  T  is  given  in  °  C  and  Psat  in  atm.  Note  that  when  X  >  1 
(i.e.  Eq.  (43))  the  water  vapour  is  in  supersaturated  conditions. 

The  single-phase  assumption  for  water  transport  in  this  model 
implies  that  water  can  exist  in  supers aturation  in  the  gas  phase 
or  if  condensed,  liquid  water  exists  as  finely  dispersed  droplets 
as  in  a  mist  flow  and  produces  negligible  effects  on  the  gas  phase 
transport.  That  is,  the  water  activity  calculated  based  on  water 
partial  pressure  is  allowed  to  be  greater  than  unity. 

The  proton  conductivity  in  the  PEM  is  dependent  upon  water 
content  X  and  temperature.  Springer  et  al.  [22]  measured  the 
proton  conductivity  of  Nation  for  a  range  of  water  activities  and 
temperatures  with  fully  hydrated  membranes.  They  found  an 
activation  energy  value  for  these  measurements.  This  activation 
energy  was  then  assumed  to  apply  for  all  values  of  X.  They 
obtained  the  following  relation  from  a  series  of  experiments 
crH+  as  a  function  of  X\ 


crH+  =  (0.00431 19A. 


0.00326)  exp 


(46) 


The  properties  reported  by  Springer  et  al.  [22]  are  used  in 
the  baseline  calculations.  The  conservation  equation  for  pro¬ 
ton  transport  in  the  PEM  is  then  given  by  Eq.  (47)  as  a  pure 
electrical  conduction  process. 


d 

dXj 


=  0 


(47) 


The  diffusion  transport  of  the  dissolved  water  in  the  mem¬ 
brane  is  related  to  its  concentration  gradient  by  an  ordinary 
diffusion  process,  i.e.  Fick’s  law  of  diffusion 


/h2Oj  =  —  7h2o  ,m 


3Ch2o 

dXj 


(48) 


where  /  H2o,m  is  the  water  diffusion  coefficient  in  the  electrolyte 
given  in  Eq.  (49) [28]: 


^H2o,m  =  1  x  10 


-10 


2.8628  In 


1818 


-  1.63795 


(49) 


The  electroosmotic  drag  is  associated  with  the  protonic 
charge.  It  has  been  shown  that  the  protonic  charge  in  the  mem¬ 
brane  has  a  drag  effect  with  water  which  is  related  to  the  polar 
structure  of  the  water  molecules  [27].  This  “adherence  effect” 
makes  the  water  molecules  stick  to  the  protons  that  are  mov¬ 
ing  from  one  charged  site  to  another,  driving  the  water  from  the 


anode  towards  the  cathode.  The  electroosmotic  drag  coefficient 
“ftd”  is  the  number  of  water  molecules  carried  per  proton  across 
the  membrane  as  electric  current  is  passed  under  conditions 
of  no  water  concentration  gradient.  Because  this  mechanism  is 
strongly  influenced  by  the  protonic  charge,  its  transport  can  be 
related  to  the  local  current  density  via  Eq.  (50) 


^H2odrag  =  nd 


(50) 


We  assume  that  electroneutrality  prevails  inside  the  PEM. 
That  is,  the  proton  concentration  in  the  polymer  is  assumed  to  be 
constant  and  equal  to  the  concentration  of  the  fixed  sulfonic  acid 
groups.  Therefore  the  proton  transport  equation  is  conserved 
inside  the  PEM.  Further  the  electroosmotic  drag  coefficient  is 
assumed  to  be  constant  and  equal  to  one  [27].  In  our  model 
electroosmosis  only  occur  where  the  local  current  density  is 
variable,  that  is  in  the  polymer  fraction  of  the  CLs,  as  indicated 
by  Eq.  (50).  This  does  not  mean  that  the  electroosmotic  drag 
has  no  effect  on  the  water  transport  through  the  membrane.  The 
electroosmotic  drag  term  enters  the  water  conservation  equation 
as  a  boundary  condition  for  the  electrolyte  fraction  at  the  CLs. 
The  total  water  flux  should  be  equal  to  the  electroosmotic  drag 
and  diffusive  flux. 

Water  absorption  and  desorption  take  place  on  the  polymer 
fraction  (25%)  inside  the  CLs.  The  rate  of  water  absorp¬ 
tion/desorption  fiH2oabs ’ *s  based  on  the  difference  between  the 
equilibrium  water  content  in  the  polymer  fraction  of  the  CLs, 
CH2o,m ,  and  the  water  content  of  the  gas  mixture  inside  the  pores 
of  the  CLs,  Cn2o,g-  The  water  content  of  either  can  be  calculated 
as 


Ch2o  = 


Pm, dry  ^ 
47m,  dry 


(51) 


where  pm,dry  and  Mm^ry  are  the  density  and  molecular  mass  of 
a  dry  membrane. 

Relation  (52)  states  that  when  the  water  content  of  the  poly¬ 
mer  fraction  is  larger  than  that  of  the  surrounding  gases,  then 
water  should  desorb  from  the  membrane.  It  also  indicates  that  if 
the  water  content  of  the  surrounding  gas  is  greater  than  that  of 
the  membrane,  then  water  will  be  absorbed  by  the  polymer. 


^H2Oabs  =  f  (CH20,m  -  CH20,g) 


(52) 


where  xj/  is  a  proportionality  constant  set  to  2.5.  Higher  or  lower 
proportionality  constants  increase  the  concentration  differences 
specially  at  high  current  densities.  The  value  of  xj/  will  mostly 
show  its  effect  under  time  dependent  calculation  (STP  or  freez¬ 
ing  conditions),  Shimpalee  et  al.  [18]  reported  values  of  ^  =1.5 
for  “start-up”  calculation  at  298  K  and  semi-dry  membranes. 

The  absorbed  and  desorbed  water  plus  the  contribution  of  the 
electroosmotic  drag  at  the  anode  (Eq.  (54))  and  cathode  (Eq. 
(53))  CLs,  are  included  into  the  water  conservation  equations 
via  source  terms 


SH20,C  =  "HzOdrag.C  +  «H2Oabs,C 


SH20,a 


+  +  «H2Oabs,a 


(53) 


(54) 
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The  first  term  on  the  right-hand- side  of  Eq.  (54)  accounts  the 
“mass  defect”  due  to  the  reactions.  Note  that  the  source  terms 
on  the  species  conservation  equations  do  not  inject  or  withdraw 
mass  but  instead  adjust  the  concentration  of  species. 

For  heat  transfer  purposes,  the  membrane  is  considered  a 
conducting  solid  slab.  The  conservation  equation  for  energy  in 
this  region  is  given  as 


d 

dxj 


,eff 


sh,  mem 


(55) 


where  Sh, mem  is  an  energy  source  term  and  k^em  is  the  effective 
heat  conductivity  of  the  membrane. 

There  is  a  heat  source  at  the  electrolyte  given  by  Eq.  (56), 
effectively  due  to  ohmic  heating 


sh,m  — 


_eff 

JH+ 


(56) 


bly. 


where  i  is  the  spatial  electric  current  value. 

2.5.  Computational  procedure 

In  most  PEM  fuel  cell  models  the  catalytic  activity  is  a 
function  of  the  concentration  of  the  reactants.  This  depen¬ 
dence  assumes  a  constant  activation  overpotential  in  the  CLs. 
This  method  makes  the  solution  process  “easy”  to  handle;  by 
assuming  a  constant  activation  overpotential,  an  average  current 
density  can  be  prescribed  before  the  calculations.  A  specific  cur¬ 
rent  density  in  the  cell  can  be  scaled  with  concentration  to  match 
a  desired  average  value. 

Bang  [2]  reported  a  more  detailed  method  regarding  the  kinet¬ 
ics  in  the  catalyst  layers.  In  this  method  the  catalytic  activity  is 
a  function  of  both  concentrations  and  spatial  variation  of  the 
overpotentials.  The  overpotentials  are  functions  of  the  concen¬ 
tration  of  species  and  overpotentials  at  their  spatial  locations. 
This  model  assumes  a  variable  overpotential  in  the  CLs,  sim¬ 
ulating  a  variable  kinetic  activity  in  these  regions.  This  way  to 
address  the  cell  overpotentials  does  not  allow  the  user  to  pre¬ 
dict  a  given  current  density  before  the  calculations.  Instead  this 
comes  as  part  of  the  solution. 

The  cell  potential  calculation  is  given  as 

^cell  =  ^rev  —  Oact  ~  Oohm  (57) 

The  overpotentials  in  Eq.  (57)  are  functions  of  their  spatial 
location.  To  be  able  to  calculate  them  at  any  given  location  in 
the  ME  A  it  was  necessary  to  define  voltage  reference  points. 
By  convention  the  land  area  of  the  anode’s  GD  was  arbitrar¬ 
ily  fixed  to  a  zero  potential  boundary  condition.  This  implies 
that  the  absolute  potential  in  the  anode’s  electrode  is  negative 
because  the  electric  charge  is  being  transferred  from  the  anode 
reaction  sites  to  the  anode’s  GD/CC  interface.  The  calculated 
cell  potential  is  given  at  the  cathode’s  GD/CC  interface.  The 
maximum  potential  is  then  given  in  the  cathode’s  CL,  as  shown 
in  Fig.  3.  Within  this  context  the  ohmic  overpotential  r]ohm  can  be 
calculated  by  the  protonic  and  electric  overpotentials  less  their 


respective  reference  potential,  as 

tarn  =  (|»7e-|  -  l^e-l)  +  (W  I  ~  I^H+l)  (58) 

The  protons  generated  in  the  anode  CL  travel  towards  the 
cathode  CL.  The  potential  losses  of  this  process  will  be  pos¬ 
itive  in  the  anode  because  of  its  positive  source  term,  su+  a, 
and  will  become  negative  as  they  travel  towards  the  cathode, 
because  of  its  negative  source  term  sH+  a.  However,  the  calcu¬ 
lations  deal  with  absolute  values  and  this  implies  that  there  will 
be  an  “assumed”  floating  zero  potential  reference  value  lying 
at  one  point  between  the  anode  and  cathode  CLs.  The  position 
of  this  reference  value  becomes  a  function  of  the  activity  of  the 
CLs.  The  overpotentials  at  each  CL  are  updated  on  an  iterative 
basis.  As  a  result  it  is  possible  to  determine  a  variable  catalytic 
activity  throughout  the  catalyst  layers.  By  separating  the  losses 
at  the  CLs  the  activation  overpotentials  were  calculated  using 

7? act,  (a, c)  =  7?tot,(a,c)  —  l^e-,(a,c)l  —  I^H+,(a,c)l  (59) 

2.6.  Solution  algorithm 

The  conservation  equations  given  in  this  paper  are  not  solved 
by  the  Star-CD  solver  in  their  differential  form,  but  by  the  finite 
volume  method,  following  a  fairly  well  known  CFD  algorithm, 
known  as  SIMPLE  [26]. 


Table  1 

Sources  and  sink  terms  of  the  conservation  equations 


Equation 

Anode  CL 

Electrolyte 

Cathode  CL 

Mass 

Sp,  a 

- 

sp,  c 

Energy 

Sh,  a 

sh,mem 

Sh,  c 

Oxygen 

- 

- 

so2 

Hydrogen 

5h2 

- 

- 

Water  vapour 

SH20,a 

~ 

5H20,c 

Potons 

sU+,a 

- 

sn+,c 

Electrons 

se~ 

~ 

se~ 
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Fig.  4.  Polarization  curves  of  the  “base  case”  and  validation  with  a  different 
numerical  model  and  experimental  data. 


Table  2 


Geometric  and  operational  parameters  used  in  the  model’s  “base  case” 


Parameter 

Value 

Unit 

Channel  length 

10 

mm 

Channel  height 

0.5 

mm 

Channel  width 

1 

mm 

GDL  thickness 

0.3 

mm 

GDL  land  area 

10 

mm2 

CL  thickness 

20 

[xm 

PEM  thickness 

50 

fxm 

Temperature 

353 

K 

Pressure 

1 

atm 

Hydrogen  stoichiometry  ratio 

2 

- 

Hydrogen  R.H.  at  inlet 

100 

% 

Air  stoichiometry  ratio 

2 

- 

Air  R.H.  at  inlet 

100 

% 

Table  3 

Base  case  membrane  electrode  assembly  parameters 

Parameter 

Symbol 

Value 

Unit 

GDL  porosity 

£gdl 

55 

% 

GDL  tortuosity 

/GDL 

1.5 

- 

GDL  heat  conductivity 

xeff 

A,mem 

0.5 

W  (rnK)"1 

GDL  electric  conductivity 

^e-,GDL 

825 

sm_1 

CL  porosity 

£cl 

13 

% 

CL  tortuosity 

/cl 

2.5 

- 

CL  heat  conductivity 

xeff 

^CL 

0.82 

W  (mK)-1 

CL  electric  conductivity 

^e-,CL 

250 

S  m_1 

Pt  mass  loading  per  unit  area  of  CL 

ra/pt 

0.6 

mg  cm-2 

Mass  of  Pt  supported  on  carbon  black 

m/7Pt 

20 

% 

Nation  volume  fraction  in  CL 

- 

25 

% 

PEM  heat  conductivity 

£eff 

^CL 

0.64 

W  (mK)-1 

Anode’s  transfer  coefficient 

Oi  a 

0.5 

- 

Cathode’s  transfer  coefficient 

Ctc 

0.5 

- 

Anode’s  exchange  current  density 

l0,a 

8.0  x  103 

Am-2 

Cathode’s  exchange  current  density 

io,c 

4.0  x  10- 

8  Am"2 

(a)  TEMPERATURE  [K] 
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Fig.  5.  Temperature  distribution  across  the  MEA  at  different  cell  potentials,  (a) 
Temperature  distribution  across  the  MEA  at  Eceii  ^  0.8  V  and  i  ~  0.01  A  cm-2 
(Tljps.eps).  (b)  Temperature  distribution  across  the  MEA  at  Eceii  &  0.35  V  and 
i  ~  0.75  A  cm-2. 


Fig.  6.  Oxygen  and  hydrogen  mass  fractions  at  i  =  0.012  Acm_2and  Eceii  = 
0.8  V.  (a)Oxygen  mass  fraction  in  the  cathode  at  0.012  Acm_2and  0.8  V.  (b) 
Hydrogen  mass  fraction  in  the  anode  at  0.012  A  cm-2  and  0.8  V. 
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The  user-defined  code  forms  part  of  this  algorithm.  It  begins 
with  a  prescribed  total  overpotential  at  the  cathode  CL  then 
the  current  density  is  first  calculated  and  stored  for  the  cath¬ 
ode  side.  Following  this,  the  anode’s  activation  overpotential  is 
tuned  in  order  to  match  the  cathode’s  total  current.  The  code 
loops  through  the  anode  catalyst  layer  a  certain  number  of  times 
per  iteration  in  order  to  approximate  the  two  total  currents.  If 
both  anode  and  cathode  currents  still  do  not  match  after  this 
number  of  loops  the  algorithm  continues  and  iterates  again.  The 
iterative  process  changes  several  variables  such  as  temperature, 
density,  diffusion  coefficients,  scalar  concentrations,  proton  con¬ 
ductivity,  overpotentials,  etc.  With  the  calculated  total  current, 
the  cell  potential  can  be  calculated  with  an  updated  value  of 
the  overpotentials.  The  procedure  continues  until  the  anode  and 
cathode  total  currents  match  (Table  1). 

The  computations  were  preformed  on  a  machine  with  a  1 .4 
GHz  processor  and  512  MB  of  memory.  To  obtain  a  single 
solution,  about  four  to  eight  thousand  iterations  were  required 
and  took  from  4  to  10  h  of  CPU  time  (the  latter  corresponds  to 
the  high  and  very  low  current  densities).  The  flowchart  diagram 
of  the  algorithm  is  shown  in  Fig.  1. 

3.  Results  and  discussions 

The  model  was  tested  under  a  range  of  operational  conditions. 
Fig.  4  shows  a  good  agreement  between  the  “base  case”  results, 


experimental  tests  performed  by  Su  et  al.  [24]  and  another  PEM 
fuel  cell  model  presented  by  Carnes  and  Djilali  [7]. 

3.1.  Base  case  model 

The  base  case  model  presented  in  this  section  is  the  foun¬ 
dation  of  the  analyses  performed  in  this  paper.  It  represents  the 
closest  attempt  by  the  authors  to  resemble  a  functional  PEM  fuel 
cell.  The  geometric  parameters  and  operation  conditions  of  the 
base  case  model  are  listed  in  Table  2. 

Table  3. 

The  temperature  field  measured  in  the  PEM  is  shown  in  Fig.  5 
for  two  different  potentials.  These  figures  show  how  the  temper¬ 
ature  can  increase  by  up  to  5  °  C  at  the  given  potentials.  In  a  PEM 
fuel  cell  the  major  heat  source  comes  from  the  entropy  change 
of  the  oxygen  reduction  reaction.  That  is  why  it  is  always  the 
top  part  (cathode)  of  Fig.  5  is  slightly  warmer  than  the  bottom 
(anode).  The  secondary  source  of  heat  in  this  type  of  cells  comes 
from  the  ohmic  resistance  of  the  PEM  to  ionic  conductivity,  this 
also  contributes  to  the  heating  of  the  cell. 

Figs.  6  and  7  show  contour  planes  used  to  illustrate  the  con¬ 
centration  of  the  reactants.  Each  plane  includes  the  channel  and 
GD  for  both  anode  and  cathode  sides  of  the  cell.  The  base  case 
model  oxygen  and  hydrogen  concentrations  are  shown  in  these 
figures  for  two  different  potentials.  The  figures  show  how,  at 
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Fig.  7.  Oxygen  and  hydrogen  mass  fractions  at  l  =  0.77  A  cm-2  and  £ceii  = 
0.37  V.  (a)  Oxygen  mass  fraction  in  the  cathode  at  0.77  A  cm-2  and  0.37  V.  (b) 
Hydrogen  mass  fraction  in  the  anode  at  0.77  A  cm-2  and  0.37  V. 
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Fig.  8.  Activation  overpotential  distributions  in  the  cathode’s  catalyst  layer  at 
different  cell  potentials,  (a)  Activation  overpotential  distribution  in  the  cathode’s 
catalyst  layer  at  0.01  A  cm_2and  0.8  V.  (b)  Activation  overpotential  distribution 
in  the  cathode’s  catalyst  layer  at  0.77  A  cm_2and  0.37  V. 
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high  current  densities,  the  reactants  are  consumed  faster  than  at 
lower  current  densities.  It  is  clear  that  at  these  current  densities 
the  most  affected  species  is  the  oxygen.  In  fact,  at  these  power 
densities  part  of  the  model’s  ME  A  is  depleted  of  oxygen  specif¬ 
ically  in  the  regions  under  the  current  collectors,  as  shown  in 
Fig.  7(a).  This  is  partly  due  to  a  limited  oxygen  diffusivity  and 
also  due  to  the  difficulty  in  removing  the  water  produced  by  the 
oxygen  reduction  reaction. 

On  the  other  hand,  it  can  be  seen  that  on  the  anode  side  of 
the  cell,  the  hydrogen  concentration  experiences  little  change. 
Hydrogen  has  a  higher  diffusivity  than  oxygen  and  this  allows  a 
better  transport  even  at  high  current  densities  and  low  porosities. 
In  the  model  presented  in  this  paper  hydrogen  posed  no  depletion 
threat.  In  addition,  the  anode  side  of  the  cell  rarely  experiences 
flooding,  therefore  its  porous  medium  is  almost  never  blocked. 

One  important  feature  of  the  model  presented  in  this  paper  is 
the  way  the  activation  overpotential  is  calculated.  As  explained 
earlier  in  Section  2.5,  the  overpotentials  of  the  cell  are  calculated 
as  a  function  of  the  species  concentration. 

The  CCL  alone,  is  represented  in  Figs.  8  and  9  and  because 
it  is  only  20  jxm  thick  it  has  been  scaled  by  several  factors  in  the 
Y  direction.  These  figures  illustrate  the  activation  overpotential 
and  the  local  current  density  distribution  in  this  region.  Fig.  8(a) 
shows  how  at  a  low  current  density  the  activation  overpotential 
is  fairly  uniform  throughout  the  CCF.  However,  as  the  current 
density  increases  the  concentration  of  the  reactants  decreases 


i  [A  cmA-2] 


Fig.  9.  Current  density  distributions  in  the  cathode’s  catalyst  layer  at  different 
cell  potentials,  (a)  Current  density  distribution  in  the  cathode’s  catalyst  layer  at 
0.01  A  cm-2  and  0.8  V.  (b)  Current  density  distribution  in  the  cathode’s  catalyst 
layer  at  0.77  A  cm-2  and  0.37  V. 


and  this  has  a  negative  effect  in  the  activation  overpotential  dis¬ 
tribution,  as  shown  in  Fig.  8(b).  Here  the  activation  overpotential 
is  not  uniform,  the  maximum  values  of  the  activation  overpo¬ 
tential  move  towards  the  lower  region  of  the  CCF  (close  to  the 
CCF/PEM  interface)  under  the  current  collector  ribs.  The  min¬ 
imum  values  always  stay  on  top  of  the  CCF  right  below  the 
cathode  channel. 

These  activation  overpotential  distributions  have  other  effects 
in  the  cell.  The  Butler- Volmer  Eqs.  (3)  and  (4)  described  the 
relation  between  the  activation  overpotential  and  the  local  cur¬ 
rent  density.  Fig.  9  shows  how  the  local  current  density  varies 
throughout  the  CCF.  At  a  low  potential  the  local  current  density 
is  evenly  distributed,  as  shown  in  Fig.  9(a).  At  high  power  densi¬ 
ties  however,  the  current  distribution  is  totally  uneven  (Fig.  9(b)). 
The  CCF  becomes  a  region  of  contrasts.  Near  the  CCF/PEM  the 
current  is  10  times  bigger  than  its  immediate  upper  region.  This 
means  that  only  the  very  bottom  layer  of  the  CCF  is  powering 
the  cell.  At  high  current  densities  the  overpotentials  near  the 
CCF/PEM  interface  are  so  big  that  almost  90%  of  the  whole 
power  density  generated  in  the  cathode  is  produced  at  this  very 
thin  layer. 

Due  to  the  importance  of  oxygen  in  the  overall  cell  perfor¬ 
mance  its  concentration  in  the  CCF  is  shown  in  Fig.  10.  Note  how 
under  the  current  collector  ribs  the  concentration  of  oxygen  is 
reduced  gradually  as  the  current  density  increases.  Once  again, 
the  same  tendency  of  the  current  density  and  activation  over- 


Fig.  10.  Oxygen  concentrations  in  the  cathode’s  catalyst  layer  at  different 
cell  potentials,  (a)  Oxygen  concentration  in  the  cathode’s  catalyst  layer  at 
0.01  A  cm_2and  0.8  V.  (b)  Oxygen  concentration  in  the  cathode’s  catalyst  layer 
at  0.77  A  cm-2and  0.37  V. 
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potential  is  reflected  in  these  figures.  At  low  current  densities 
the  availability  of  oxygen  is  good  enough  to  fulfill  the  reactant 
requirements  (Fig.  10(a)).  However,  at  high  current  densities  the 
oxygen  diffusivity  becomes  a  problem,  as  it  cannot  provide  the 
rate  needed  to  drive  the  reaction  further.  At  this  current  densities, 
starvation  regions  are  clearly  present  under  the  current  collector 
ribs,  Fig.  10(b). 

As  explained  earlier  in  Section  2.5,  the  activation  overpoten¬ 
tial  is  a  sum  of  other  overpotentials  of  a  different  nature.  The 
following  figures  show  the  distribution  of  these  losses  in  the  cell. 
Ohmic  overpotentials:  ionic  and  electronic,  are  always  present 
and  play  active  roles  in  the  calculation  of  the  activation  overpo¬ 
tential.  Due  to  the  intrinsic  electric  resistance  of  the  materials  in 
the  cell  the  losses  of  these  two  overpotentials  become  greater  as 
the  current  density  increases. 

Fig.  1 1  shows  the  protonic  overpotential  distribution  at  two 
different  current  densities.  The  loss  is  the  most  severe  out  of 
the  two  ohmic  overpotentials.  In  this  model  it  is  a  function  of 
the  membrane  water  content  and  temperature.  The  proton  and 
water  distribution  in  the  membrane  are  related  by  the  electroos- 
motic  drag  effect.  A  special  procedure  to  model  this  relationship 
was  presented  in  Section  2.5.  This  method  requires  a  reference 
potential  (0  V)  located  in  between  the  two  CLs.  However,  this 
reference  potential  has  a  variable  position;  Fig.  1 1(a)  shows  how 
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Fig.  11.  Protonic  overpotential  distributions  at  different  cell  potentials,  (a) 
Protonic  overpotential  distribution  at  0.01  Acm_2and  0.8  V.  (b)  Protonic  over¬ 
potential  distribution  at  r]H+c  at  0-77  A  cm-2  and  0.37  V. 


Fig.  12.  Electric  overpotential  distributions  at  different  cell  potentials,  (a) 
Electric  overpotential  distribution  at  0.77  Acm_2and  0.38  V.  (b)  Electric  over¬ 
potential  distribution  at  0.77  A  cm_2and  0.37  V. 


at  a  low  current  is  hard  to  distinguish  were  the  reference  value 
lies.  At  this  power  density  the  ohmic  losses  are  small  and  the 
biggest  losses  are  located  at  the  anode  CL.  Fig.  11(b)  shows 
how  the  reference  value  lies  clearly  in  between  the  two  CLs 
adopting  a  curved  shape  under  the  current  collector  ribs.  The 
biggest  losses  at  this  power  density  are  clearly  located  close  to 
the  cathode  CL. 

Fig.  12  shows  the  distribution  of  the  electric  overpotentials  in 
the  GDs.  Similar  to  the  protonic  losses,  the  electric  overpotential 
also  uses  reference  voltages.  In  this  case  these  reference  poten¬ 
tials  are  stationary.  These  are  set  as  boundary  conditions  at  the 
interfaces  between  the  GDs  and  the  current  collector  ribs.  This 
condition  is  reflected  in  the  shape  of  the  contours  of  Fig.  12.  The 
electric  overpotential  is  highest  at  the  region  under  the  channels 
and  closer  to  the  inlet  boundaries.  Once  again,  it  can  be  seen  that 
the  higher  the  current  density  the  higher  the  overpotential. 

The  protonic  and  electric  overpotential  are  two  conduction 
processes  restricted  by  the  resistance  of  the  materials.  The 
protonic  diffusion  coefficient  is  calculated  as  a  function  of  tem¬ 
perature  and  water  content  of  the  membrane  (see  Section  2.4.4). 


4.  Conclusions 


A  comprehensive  three-dimensional  model  of  a  PEM  fuel 
cell  has  been  developed.  The  model  accounts  for  all  major  trans¬ 
port  phenomena  in  the  cell  with  a  single-phase  assumption.  The 
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transport  equations  and  boundary  conditions  of  each  of  the  cell’s 
regions  were  introduced  along  with  some  special  boundary  con¬ 
ditions  such  as  the  inlet  velocity  and  the  CL  and  PEM  interfaces. 
A  thorough  explanation  of  the  processes  occurring  inside  the 
electrolyte  was  given  in  Section  2.4.4.  It  began  with  the  equi¬ 
librium  water  content  and  continued  with  water  diffusion  and 
electroosmotic  drag  and  their  impact  upon  proton  and  water 
conductivity. 

We  produced  a  good  agreement  of  the  model’s  performance 
with  some  other  experimental  and  numerical  models.  Some  illus¬ 
trations  of  the  model  operating  at  different  conditions  showed 
how  the  reactants  distribute  throughout  the  cell.  The  temperature 
distribution  in  the  PEM  at  different  cell  potentials  showed  how 
the  temperature  gradient  could  hinder  the  material  integrity  of 
the  membrane.  It  was  shown  how  the  current  density  distributes 
across  the  cathode  CL  (Fig.  9).  To  our  knowledge  these  current 
density  distributions  have  long  been  measured  in  experiments 
but  never  published  in  a  computational  model.  Similarly  the  dis¬ 
tribution  of  the  overpotentials  were  shown  in  the  gas  diffusers, 
catalyst  layers  and  in  the  electrolyte. 
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